package scu.maqiang.numeric;

import java.util.Arrays;
import java.util.*;

/**
 * 
 * @author Administrator
 *
 */
public class EigenDSolver extends Solver {
	public EigenDSolver() {

	}

	public EigenDSolver(double eps, int maxIteration) {
		this.eps = eps;
		this.maxIteration = maxIteration;
	}

	/**
	 * 
	 * @param A      稠密矩阵A
	 * @return 特征值与对应的特征向量
	 */
	public Pair<Double, double[]> powerit(double[][] A) {
		int n = A.length;
		ParamCheck.checkEqual(n, A[0].length);
		double lambda = 0.0;
		double[] x = new double[n];
		Random random = new Random(System.currentTimeMillis());
		Arrays.setAll(x, i -> random.nextDouble());
		// System.out.println(Arrays.toString(x));
		double lambda0 = Double.MAX_VALUE;
		double L2Norm = 0.0;
		double[] Ax = x.clone();
		for (int i = 0; i < maxIteration; i++) {
			L2Norm = MVO.L2Norm(x);
			for (int j = 0; j < n; j++) {
				x[j] = x[j] / L2Norm;
			}
			MVO.matmul(A, x, Ax);
			lambda = MVO.dot_product(Ax, x);
			if (Math.abs(lambda - lambda0) < eps) {
				break;
			}
			lambda0 = lambda;
			for (int j = 0; j < n; j++) {
				x[j] = Ax[j];
			}
		}
		return new Pair<Double, double[]>(lambda, x);
	}

	public void inversePowerTran(double[][] A, double[][] M, double[] x0, double alpha, int num, double[] D,
			double[][] V) {
		int i, j, k;
		int n = A.length;
		double[] x = new double[n];
		double[] xx = new double[n];
		double[] y = new double[n];
		double[] summ = new double[n];
		double r = 0.0, rr = 0.0;
		double KK = 0.0, MM = 0.0;
		double[][] L = new double[n][n];
		for (i = 0; i < n; i++) {
			for (j = 0; j < n; j++) {
				L[i][j] = A[i][j] + alpha * M[i][j];
			}
		}
		DirectDSolver solver = new DirectDSolver(L);
		solver.PLU();
		for (i = 0; i < num; i++) {
			for (int kk = 0; kk < n; kk++) {
				x[kk] = x0[kk];
			}
			k = 1;
			while (k < 100) {
				Arrays.fill(summ, 0.0);
				for (j = 0; j < i; j++) {
					double[] Mx = MVO.matmul(M, x);
					double temp = MVO.dot_product(V[j], Mx);
					for (int jj = 0; jj < n; jj++) {
						summ[jj] += V[j][jj] * temp;
					}
				}
				for (int ii = 0; ii < n; ii++) {
					x[ii] -= summ[ii];
				}
				double[] Mx = MVO.matmul(M, x);
				double xMx = Math.sqrt(MVO.dot_product(x, Mx));
				for (int ii = 0; ii < n; ii++) {
					x[ii] /= xMx;
				}
				y = MVO.matmul(M, x);
				double[] yy = y.clone();
				solver.PLUSolve(y, xx);
				//System.out.println("Find\t" + (i + 1) + "\t of\t" + num + "\t eigenvalues!");
				KK = MVO.dot_product(xx, yy);
				double[] Mxx = MVO.matmul(M, xx);
				MM = MVO.dot_product(xx, Mxx);
				rr = KK / MM;
				if (k > 1 && Math.abs((rr - r) / rr) < this.eps) {
					break;
				}
				k++;
				r = rr;
				for (int kk = 0; kk < n; kk++) {
					x[kk] = xx[kk];
				}

			}
			//System.out.println(i + "  of  " + num + "  finished!");
			// write(*, *) rr - alpha
			D[i] = rr - alpha;
			for (int kk = 0; kk < n; kk++) {
				V[i][kk] = xx[kk] / Math.sqrt(MM);
			}
		}
	}
}
